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Abstract. 

Characteristic methods show excellent promise in the evolution of single 
black hole spacetimes. The effective coupling with matter fields may help the 
numerical exploration of important astrophysical systems such as neutron star black 
hole binaries. To this end we investigate formalisms for numerical relativistic 
hydrodynamics which can be adaptable to null (characteristic) foliations of the 
spacetime. The feasibility of the procedure is demonstrated with one-dimensional 
' results on the evolution of self-gravitating matter accreting onto a dynamical black 

hole. 
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Introduction 

The so-called characteristic formulation of the field equations of general relativity has 
proven recently to be well suited for numerical evolutions of single black hole (vacuum) 
spacetimes |l| . A three-dimensional characteristic code developed by the Binary Black 
Hole Grand Challenge Alliance in the US has achieved, for the first time, accurate 
and long-term stable evolutions of such spacetimes, including distorted, moving and 
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spinning single black holes, with evolution times up to 60, OOOM Such evolutions 
are ultimately possible due to the excellent gauge control properties furnished by null 
foliations, gauge control being a major open issue in the more usual 3+1 formulation. 

The incorporation of matter terms into such a computational framework, in 
the form of some suitable stress-energy tensor, permits the study of interesting 
astrophysical scenarios with optimal computational efficiency. Such an approach 
allows to investigate the non-linear evolution of self-gravitating matter around dynamic 
black holes. Some interesting astrophysical studies to be performed include, e.g., 
accretion disks or the non-linear and fully relativistic evolution of a compact binary 
system formed by a black hole and a neutron star. This binary system is one of 
the anticipated sources of gravitational radiation. Its detailed numerical study can 
provide valuable information on the waveform signals and a means to decipher the 
nuclear equation of state. 

In this report we focus on the formulation of the hydrodynamic equations on null 
foliations of the spacetime. The reader must be aware that most of the previous 
analytic and numerical work in general relativistic hydrodynamics (GRH hereafter) 
has been done using the "standard" 3+1 formulation, i.e., using a spacetime foliation 
of spacelike hypersurfaces (see, e.g., || and references therein). Additionally, the 
small amount of existing work dealing with the numerical integration of the GRH 
equations on null hypersurfaces can be characterized by: a) a direct way of writing 
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the hydrodynamic equations with no attention to important mathematical (and 
numerical) details (as, e.g., hyperbolicity) and b) the use of straightforward finite- 
differences numerical schemes However, on the basis of the experience 
gained in the last few years in the integration of the hydrodynamic equations in the 
spacelike case, we believe that the procedure previously adopted for the null case is 
likely to be inadequate when dealing with non-trivial simulations. Most notably we 
refer to realistic astrophysical scenarios containing high speed (ultrarclativistic) flows 
and steep gradients or even shock waves. 

As a way to avoid these problems we propose the use of the so-called high- 
resolution shock- capturing schemes (HRSC in the following) to integrate the GRH 
equations on light cones. HRSC schemes are commonly used in numerical integrations 
of the Newtonian inviscid Euler equations. In recent years they have been 
succesfully extended to relativistic hydrodynamics, both in the special and general 
cases R,fl)H,ftl|10|- Mathematically, HRSC schemes rely on the hyperbolic 
character of the hydrodynamic equations. The knowledge of the characteristic fields 
(eigenvalues) of the system, together with the corresponding eigenvectors, allows for 
accurate integrations, by means of either exact or approximate Riemann solvers, along 
the fluid characteristics. 

Among the interesting properties of HRSC schemes we stress the important fact 
that they are written in conservation form. Hence, all physically conserved quantities 
of a given partial differential equation will also be conserved in its finite-differenced 
version. More precisely, the time update of the following hyperbolic one-dimensional 
partial differential system of equations 

du di . 

is done according to the following algorithm: 

< +1 = < - ^6+1/2 - fi-1/2) + Ais,. (2) 



In these expressions f and s are vector-valued functions, the fluxes and sources, and 
they both depend on the state vector of the system u. Index i labels the grid location 
and n the time. Time and space discretization intervals are indicated by At and Ax, 
respectively. The "hat" in the fluxes is used to denote the so-called numerical fluxes 
which, in a HRSC scheme, are computed according to some generic flux-formula. It 
adopts the following functional form: 

l±i = \ (f(uf ± .) + f(uf ±| ) - £ I K I AZ a r a ] . (3) 

Notice that the numerical flux is computed at cell interfaces (i ± 1/2). Indices L and 
R indicate the left and right sides of a given interface. The sum extends to p, the 
total number of equations. Finally, quantities A, Aw and r denote the eigenvalues, 
the jump of the characteristic variables and the eigenvectors, respectively, computed 
at the cell interfaces according to some suitable average of the state vector variables. 
Further technical information can be found in, e.g. || and references therein. 

HRSC schemes are also known for giving stable and sharp discrete shock profiles. 
They have also a high order of accuracy, typically second order or more, in smooth 
parts of the solution. 



3 



2. A covariant form of the general relativistic hydrodynamic equations 

The conservation equations in covariant form are 



— y /=gr^ = y/=gr^ (5) 
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with the matter current and stress energy tensor for a perfect fluid given by J M = pu^, 
T^ v = phu^u w + pg^ v , where p is the density, p the pressure and u^ 1 the fluid four 
velocity obeying the normalization condition u^u^ = —1. Additionally, h = 1 + e+p/p 
is the relativistic specific enthalpy, with e being the specific internal energy 

The construction of an initial value problem requires the introduction of a 
coordinate chart (x°,x l ), where the scalar x° represents a foliation of the spacetime 
with hypersurfaces (coordinatized by (x , x 2 , x 3 )) on which an appropriate initial data 
set is prescribed. Upon introducing the coordinate chart, we define the coordinate 
components of the four- velocity u M = (u°,u l ). The velocity components u l , together 
with the rest-frame density and internal energy, p and e, provide a unique description 
of the state of the fluid and are called the primitive variables. They consitute a vector 
in a five dimensional space w A — (p, u l , e). The index A is taken to run from zero to 
four, coinciding for the values (1,2,3) with the coordinate index i. 

We define the initial value problem in terms of another vector in the same fluid 
state space, namely the conserved variables U s , individually denoted (D,S l ,E), as 
follows: 

D = U° = J° = pu" (6) 
S* = LP = T w = phu°u l + pg Qi (7) 
E = U 4 = T m = phu°u° + pg 00 (8) 

With these definitions the equations take the standard conservation law form 

d x o(V=gU A ) + d xJ (V=gF jA ) = s A (9) 

where S' 4 (w s ) are source terms that depend on metric derivatives and the 
(undifferentiated) stress energy tensor, y/—g is the volume element associated with 
the four-metric and we define the flux vectors F- jA 

F 3 ° = J j = pv? (10) 
F ji = T ji = phu i u j + pg ij (11) 
F ji _ T j0 = phu o u j + pg 0j (12 ) 

The source terms take the form 

S° = (13) 

S l =-^r; A T^ A (14) 

S 4 = - ^—gTl x T^ (15) 

For a complete description of the local characteristic structure of the previous 
equations, which is the basic ingredient to make use of HRSC schemes, we refer the 
interested reader to . We end this section by stressing that the previous formulation 
of the GRH equations is fully covariant and, hence, both applicable to spacelike or 
null spacetime foliations. The results we show next are specialized to null foliations. 
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3. Results 

3.1. The Riemann "problem on a null surface 

This section is devoted to demonstrate the functionality of our approach in a standard 
hydrodynamical testbed computation, the Riemann problem. This is one of the 
standard tests to calibrate numerical schemes in fluid dynamics. It provides a strong 
check of our procedure in the limiting case of Minkowski spacctime. Initial data are 
given on the light cone and the evolution is compared to the exact solution. The initial 
setup consists of a zero velocity fluid having two different thermodynamical states on 
either side of an interface. In our setup the initial state is specified by p r = 13.3, 
p r = 10 at the right side of the interface and pi = 0.66 • 10~ 6 , pi = 1 at the left side. 
We use a perfect fluid equation of state, p — (T — l)pe with T = 5/3. When the 
interface is removed, the fluid evolves in such a way that four constant states occur. 
Each state is separated by one of three elementary waves: a shock wave, a contact 
discontinuity and a rarefaction wave. 

In Fig. |l] we plot the results of the simulation. The left pannels show the whole 
domain, with the x-coordinate ranging from to 1000 whereas the right pannels show 
a closed-up view of the most interesting region (x € [180,820]). We use a numerical 
grid of 2000 zones and, hence, a rather coarse spatial resolution (Ax — 0.5). The 
initial data are evolved up to a final time v = 120. From top to bottom Fig. [l] displays 
the internal energy (e), the velocity (u x ) and the density (p). The thick dotted line 
represents the initial discontinuity. The solid line shows the exact solution after a 
retarded time v — 120 and the thin dotted line indicates the numerical solution at this 
same time. The agreement between the two solutions is remarkable. The flow solution 
is characterized by an (ingoing) shock wave (moving to the left) and an (outgoing) 
rarefaction wave (moving to the right). We note that those features of the solution 
moving towards the left are more developed than those moving towards the right. 
The major numerical deviations from the exact solution appear at the post contact 
discontinuity region, with the density and the internal energy being slightly under and 
over estimated there. 



3.2. Spherical accretion onto a black hole 

We further test our procedure by increasing the complexity of the numerical simulation 
with the inclusion of gravity. To this end, we consider the problem of spherical 
accretion of a perfect fluid onto a black hole. 

We consider the general spherically symmetric spacetime with perfect fluid matter, 
following the formalism of Tamburino-Winicour |fl3f . The use of a spacetime foliation 
v by null hypersurfaces (V M wV M f = 0), together with the assumption of spherical 
symmetry implies that the non-trivial Einstein equations split in the following group 
of equations 

G vr — nT vr , (16) 
(jff — KiT^f , (17) 
G,„,|r = K,T vv \r, (18) 



where the last equation is a conservation condition to be imposed on the boundary T 
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Figure 1. The shock tube problem. Exact versus numerical solution at a retarded 
time v = 120 for the internal energy (top), velocity (middle) and density (bottom). 



of the integration domain. Adopting the Bondi-Sachs form of the metric element, 
p 2 PV 

ds 2 = -dv 2 + 2e 2(3 dvdr + r 2 {d9 2 + sintf 2 ^ 2 ), (19) 

r 

the geometry is completely described by the functions /3(v,r) and V(v,r). 

The first two equations above, Eqs. (16) and (17), are hypersurface equations, to 
be integrated away from T along the null cone. The (3 hypersurface equation is 

P, r = 2irrT rr (20) 
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and the V hypersurface equation, 

V r = e 2 ?(l - Anr 2 (g AB T AB - T)). (21) 

Boundary conditions /3(n)r, V{v)t for the radial integrations are provided by the 
third equation, after adopting a suitable gauge condition. This condition fixes the only 
remaining freedom in the coordinate system, namely the rate of flow of coordinate 
time at the world-tube. We refer the reader to |L2j for a discussion on suitable gauge 
conditions. 

The initial data consist of the fluid variables (p, e, u r ) on the initial slice vq, together 
with the values of /3(t>o)r & n d V(i>o)r- 

In the limit of a test fluid, i.e., a perfect fluid with sufficiently low density accreting 
spherically onto a static black hole, the numerical evolution can be compared to an 
exact solution. This comparison was performed and the agreement was found to be 
excellent. We refer the interested reader to [|l2| (see also ITlj]) for further information. 
We focus here on the case of spherical accretion onto a dynamic black hole. 
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Figure 2. Spherical accretion of a self-gravitating perfect fluid: evolution of the 
Riemann scalar curvature. Once the shell is swallowed by the central object the 
solution is dominated by the curvature of the final black hole. 

As initial data we consider the exact, stationary solution corresponding to the 
test fluid case and, on top of this background solution, we build a self-gravitating 
spherically symmetric shell whose density is parametrized by a Gaussian profile with 
suitable width and amplitude. The results of the simulation are plotted in Figs. || 
and ||. We find that the shell is radially advected (accreted) towards the hole in the 
first 10M — 15M. Once the bulk of the accretion process ends, we are left with a quasi- 
stationary background solution. This is clear in Fig. 0, where we plot the evolution of 
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the logarithm of the Riemann scalar curvature, as a function of the ingoing Eddington- 
Finkelstein radial coordinate, r. This variable encodes a large amount of relevant 
metric information. We note that the initial shell has a non-negligible curvature. 
After the shell is accreted by the hole the curvature profile monotonically increases 
towards the central singularity. The accretion of the shell originates a rapid increase 
on the mass of the apparent horizon of the black hole. This is depicted in Fig. The 
horizon almost doubles its size during the first 10M — 15M (this is enlarged in the 
insert of Fig. |3|). Once the main accretion process has finished, the mass of the horizon 
slowly increases, in a quasi-steady manner, whose rate depends on the mass accretion 
rate imposed at the world-tube, T, of the integration domain. We emphasize the fact 
that this numerical solution can be evolved as far as desired into the future. It has 
proven to remain remarkably stable and free of any kind of numerical pathologies, 
including secular instabilities or exponentially growing modes. 




Figure 3. Spherical accretion of a self-gravitating perfect fluid: evolution of the 
black hole apparent horizon mass. 
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